clear all 
set maxvar 20000

use "data\out\elbileiere2019_maalepunkt2019_wattmete2018.dta", clear 
keep Navn Merke* Modell* Måle* Modår* mergedon 
keep if inlist(mergedon, "Navn", "name1Adresse1")
drop mergedon 
g n = _n 
reshape long Merke Modell Målepunktnr Modår, i(n)
drop n _j 
drop if Målepunktnr == . 
destring Modår, g(year)
merge m:1 Merke Modell Modår using "data\in\elbiltyper_Modell_Modår_Antall_Batterikapasitet_kWh.dta", keepusing(Batterikapasitet) nogen 
replace Batterikapasitet = 16 if Batterikapasitet == . // assume missing is small battery 
su Batterikapasitet, detail 
g byte abovemedianbattery = (Batterikapasitet > `r(p50)')
g byte EVlarge = abovemedianbattery 
drop year 
tab Merke EVlarge, m 
keep Målepunktnr EVlarge Batterikapasitet
rename Målepunktnr maalepktnr
*gsort maalepktnr -EVlarge 
gsort maalepktnr -Batterikapasitet // keep larges EV battery in household 
keep if maalepktnr != maalepktnr[_n-1]
tempfile EVlarge
save `EVlarge', replace 

use "data\out\estsample.dta", clear
merge m:1 maalepktnr using `EVlarge', keep(master match) nogen 
tab  EVlarge elcar, m 
replace EVlarge = 0 if EVlarge == .
replace EVlarge = 0 if elcar == 0  

drop TNPE shoulder TsE NonPeakshoulder Effektpris-Effektprisshoulderelcar 
g byte PeakxEvent = Peak*EventDay

* day before 
g byte EventDay_pre = 0 
g byte EventDay_pre2 = 0 
g byte EventDay_pre3 = 0 
foreach d in 21893 21902 21937 21944 21958 21971 21979 {
  replace EventDay_pre = 1 if dato == `d' - 1 
  replace EventDay_pre2 = 1 if dato == `d' - 2 
  replace EventDay_pre3 = 1 if dato == `d' - 3 
}
forval h = 0/23 {
g byte Tpre`h'E = Treat * (hour == `h') * EventDay_pre * (1 - elcar)
label variable Tpre`h'E "`h'"
g byte Tpre`h'Eelcar = Treat * (hour == `h') * EventDay_pre * elcar
label variable Tpre`h'Eelcar "`h'"
}
* treatment day 
forval h = 0/23 {
g byte T`h'E = Treat * (hour == `h') * EventDay * (1 - elcar)
label variable T`h'E "`h'"
g byte T`h'Eelcar = Treat * (hour == `h') * EventDay * elcar
label variable T`h'Eelcar "`h'"
}
* day after 
g byte byte EventDay_post = 0 
foreach d in 21893 21902 21937 21944 21958 21971 21979 {
  replace EventDay_post = 1 if dato == `d' + 1 
}
forval h = 0/23 {
g byte Tpost`h'E = Treat * (hour == `h') * EventDay_post * (1 - elcar)
label variable Tpost`h'E "`h'"
g byte Tpost`h'Eelcar = Treat * (hour == `h') * EventDay_post * elcar
label variable Tpost`h'Eelcar "`h'"
}

forval h = 0(2)22 {
label variable Tpre`h'E "`h'"
label variable Tpre`h'Eelcar "`h'"
label variable T`h'E "`h'"
label variable T`h'Eelcar "`h'"
label variable Tpost`h'E "`h'"
label variable Tpost`h'Eelcar "`h'"
}
forval h = 1(2)23 {
label variable Tpre`h'E ""
label variable Tpre`h'Eelcar ""
label variable T`h'E ""
label variable T`h'Eelcar ""
label variable Tpost`h'E ""
label variable Tpost`h'Eelcar ""
}
/*
foreach y in TPnextdays TNPnextdays  {
g `y'elcar = `y'*elcar 
}
*/

keep lnforbruk T* *elcar hour temp* maalepktnr dato EventDay* EVlarge
drop TPE 
g byte weekday = dow(dato)
g byte Peak = inrange(hour,4,10)
g byte PeakxEvent = Peak*EventDay

areg lnforbruk i.dato Tpre* T0E T1E T2E T3E T4E T5E T6E T7E T8E T9E T10E T11E T12E T13E T14E T15E T16E T17E T18E T19E T20E T21E T22E T23E T?Eelcar T??Eelcar ///
  Tpost* ///
    hour hour#elcar PeakxEvent PeakxEvent#elcar /// 
  c.temp##c.temp##c.temp temp_*  c.temp#elcar ///
  if elcar == 0 | (elcar == 1 & EVlarge == 1) /// 
  /// if dato <= date("2019_12_10","YMD") + 2 ///
  , a(maalepktnr) vce(cl maalepktnr) 
est store ev2_3large 
estadd ysumm  

estout ev2_3large using "data\out\figure_C6_EVlarge.txt", replace type ///
  keep(T* PeakxEvent* *hour *elcar) ///
  cells(b( fmt(3)) se(par fmt(3)) _star) indicate(temp *dato*) mlabels(1 2 3 4 5 6 ) collabels(none) legend ///
  varlabels(TPE TP TNPE TNP) starlevels(* 0.10 ** 0.05 *** 0.01) ///
  stats(ymean r2 N, fmt(3 3 0 0)) //style(tex)
  


areg lnforbruk i.dato Tpre* T0E T1E T2E T3E T4E T5E T6E T7E T8E T9E T10E T11E T12E T13E T14E T15E T16E T17E T18E T19E T20E T21E T22E T23E T?Eelcar T??Eelcar ///
  Tpost* ///
    hour hour#elcar PeakxEvent PeakxEvent#elcar /// 
  c.temp##c.temp##c.temp temp_*  c.temp#elcar ///
  if elcar == 0 | (elcar == 1 & EVlarge == 0) /// 
  /// if dato <= date("2019_12_10","YMD") + 2 ///
  , a(maalepktnr) vce(cl maalepktnr) 
est store ev2_3small 
estadd ysumm  

estout ev2_3small using "data\out\figure_C6_EVsmall.txt", replace type ///
  keep(T* PeakxEvent*  *hour *elcar) ///
  cells(b( fmt(3)) se(par fmt(3)) _star) indicate(temp *dato*) mlabels(1 2 3 4 5 6 ) collabels(none) legend ///
  varlabels(TPE TP TNPE TNP) starlevels(* 0.10 ** 0.05 *** 0.01) ///
  stats(ymean r2 N, fmt(3 3 0 0)) //style(tex)

  
  
* Figures

graph set window fontface "Times New Roman"

* Below median EV battery size 
import delimited "data\out\figure_C6_EVsmall.txt", clear

g sample = "No and small EV"
drop if _n == 1 
rename (v1 v2) (variable b_se_star)
replace variable = variable[_n-1] if variable == ""
drop if strpos(b_se_star, "*") > 0 | b_se_star == "" 
keep if strpos(variable, "T") > 0 
g b = b_se_star[_n-1] if variable == variable[_n-1]
g se = b_se_star 
replace se = subinstr(se, "(", "", .)
replace se = subinstr(se, ")", "", .) 
keep if b != ""
drop b_se_star 
destring b se, replace 
g byte elcar = (strpos(variable, "elcar") > 0)
g day = "pre" if strpos(variable, "pre") > 0 
replace day = "post" if strpos(variable, "post") > 0 
replace day = "event" if day == "" 
tempfile small 
save `small', replace 

import delimited "data\out\figure_C6_EVlarge.txt", clear

g sample = "No and large EV"
drop if _n == 1 
rename (v1 v2) (variable b_se_star)
replace variable = variable[_n-1] if variable == ""
drop if strpos(b_se_star, "*") > 0 | b_se_star == "" 
keep if strpos(variable, "T") > 0 
g b = b_se_star[_n-1] if variable == variable[_n-1]
g se = b_se_star 
replace se = subinstr(se, "(", "", .)
replace se = subinstr(se, ")", "", .) 
keep if b != ""
drop b_se_star 
destring b se, replace 
g byte elcar = (strpos(variable, "elcar") > 0)
g day = "pre" if strpos(variable, "pre") > 0 
replace day = "post" if strpos(variable, "post") > 0 
replace day = "event" if day == "" 

append using `small'

destring variable, i(T pre elcar elcar post E) g(hour)

foreach x in b {
	gen `x'_lo = `x'-1.96*se
	gen `x'_hi = `x'+1.96*se
}
*egen yhi = max(b_hi) 
*egen ylo = min(b_lo)
g yhi = 0.2 
g ylo = -0.4 
g hour1 = hour - 0.2 
g hour2 = hour 
g hour3 = hour + 0.2 

two (line b hour1 if sample == "No and small EV" & elcar == 0, c(L) lc(black)  lp(l)) ///
  (rcap b_lo b_hi hour1 if sample == "No and small EV" & elcar == 0, lc(black) lw(vthin)) ///
  (sc b hour if sample == "No and large EV" & elcar == 1, m(diamond)  msize(vsmall) mcolor(forest_green) c(L) lc(forest_green) lp(-) ) ///
  (rcap b_lo b_hi hour if sample == "No and large EV" & elcar == 1,  lc(forest_green) lw(vthin)) ///
  (sc b hour3 if sample == "No and small EV" & elcar == 1, m(X) msize(small) mcolor(green%50) c(L) lc(green%50) lp(-) ) ///
  (rcap b_lo b_hi hour3 if sample == "No and small EV" & elcar == 1,  lc(green) lw(vthin)) ///
  (rarea yhi ylo hour1 if inrange(hour,16,22), col(gs11%20)) /// 
  if day == "event" /// 
  , xlab(0(2)23, glcolor(gs13)) ylab(, glcolor(gs13)) ytitle("") xtitle("") ///xline(16) xline(21) 
  legend(order(1 "No EV" 3 "Large EV" 5 "Small EV") col(3) pos(6)) /// 
  title("a) CPP day") /// 
  name(a, replace)

two (line b hour1 if sample == "No and small EV" & elcar == 0, c(L) lc(black)  lp(l)) ///
  (rcap b_lo b_hi hour1 if sample == "No and small EV" & elcar == 0, lc(black) lw(vthin)) ///
  (sc b hour if sample == "No and large EV" & elcar == 1, m(diamond)  msize(vsmall) mcolor(forest_green) c(L) lc(forest_green) lp(-) ) ///
  (rcap b_lo b_hi hour if sample == "No and large EV" & elcar == 1,  lc(forest_green) lw(vthin)) ///
  (sc b hour3 if sample == "No and small EV" & elcar == 1, m(X) msize(small) mcolor(green%50) c(L) lc(green%50) lp(-) ) ///
  (rcap b_lo b_hi hour3 if sample == "No and small EV" & elcar == 1,  lc(green) lw(vthin)) ///
  (rarea yhi ylo hour1 if inrange(hour,16,22), col(gs11%20)) /// 
  if day == "post" /// 
  , xlab(0(2)23, glcolor(gs13)) ylab(, glcolor(gs13)) ytitle("") xtitle("") ///xline(16) xline(21) 
  legend(order(1 "No EV" 3 "Large EV" 5 "Small EV") col(2) pos(6)) /// 
  title("b) Day after CPP day") /// 
  name(b, replace)

grc1leg2 a b, ycommon xcommon scale(1.5) name(figure_C6, replace) ysize(3)
graph export "output\figure_C6.pdf", as(pdf) name("figure_C6") replace 
